/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */

package MF2DSolver;

import HashMatrix2.HashMatrix2;
import HashMatrix2.LUFactorization;
import Integration.GaussQuadrature;
import Integration.IFunction1D;
import Integration.IFunction2D;
import MathLib.HashMatrix.Gauss;
import MathLib.HashMatrix.HashMatrix;
import MathLib.triangulation.SaveStructures.Node;
import MeshlessShapeFunctions.JGaussian;
import MeshlessShapeFunctions.JMeshlessInterpolationFunction;
import MeshlessShapeFunctions.JMeshlessInterpolationFunctionLiu;
import MeshlessShapeFunctions.Particle;
import java.util.ArrayList;
import java.util.Iterator;
import java.util.logging.Level;
import java.util.logging.Logger;
import syntaxanalyser.PostfiksXY;

/**
 *
 * @author mark_o
 */
public class MF2DSolver {
    private static final double eps = 1E-9;
	protected ArrayList<Particle> NodeList;
    protected ArrayList<JMeshlessInterpolationFunctionLiu> ListOfShapeFunction;
    protected HashMatrix2 A;
    protected double[] b;

    String pString="1";
    String wString="0";
    String qString="0";
    String fString="100";

        IFunction2D p;
        IFunction2D w;
        IFunction2D q;
        IFunction2D f;


    double left;
    double right;
    double top;
    double bottom;

    boolean DirihetAx;
    double AxVal;
    boolean DirihetBx;
    double BxVal;
    boolean DirihetAy;
    double AyVal;
    boolean DirihetBy;
    double ByVal;
    public int NNN;


    public void generate(){

        p = new IFunction2D() {

            public double calculate(double x,double y) {
                try {
                    return PostfiksXY.Calculate(pString, x,y);
                } catch (Exception ex) {
                    Logger.getLogger(MF2DSolver.class.getName()).log(Level.SEVERE, null, ex);
                }
                return 0;
            }

            public double calculateDerivatyX(double x, double y) {
                throw new UnsupportedOperationException("Not supported yet.");
            }

            public double calculateDerivatyY(double x, double y) {
                throw new UnsupportedOperationException("Not supported yet.");
            }
        };
        w = new IFunction2D() {

            public double calculate(double x,double y) {
                try {
                    return PostfiksXY.Calculate(wString, x,y);
                } catch (Exception ex) {
                    Logger.getLogger(MF2DSolver.class.getName()).log(Level.SEVERE, null, ex);
                }
                return 0;
            }

            public double calculateDerivatyX(double x, double y) {
                throw new UnsupportedOperationException("Not supported yet.");
            }

            public double calculateDerivatyY(double x, double y) {
                throw new UnsupportedOperationException("Not supported yet.");
            }
        };
        q = new IFunction2D() {

            public double calculate(double x,double y) {
                try {
                    return PostfiksXY.Calculate(qString, x,y);
                } catch (Exception ex) {
                    Logger.getLogger(MF2DSolver.class.getName()).log(Level.SEVERE, null, ex);
                }
                return 0;
            }

            public double calculateDerivatyX(double x, double y) {
                throw new UnsupportedOperationException("Not supported yet.");
            }

            public double calculateDerivatyY(double x, double y) {
                throw new UnsupportedOperationException("Not supported yet.");
            }

        };
        f = new IFunction2D() {

            public double calculate(double x,double y) {
                try {
                    return PostfiksXY.Calculate(fString, x,y);
                } catch (Exception ex) {
                    Logger.getLogger(MF2DSolver.class.getName()).log(Level.SEVERE, null, ex);
                }
                return 0;
            }

            public double calculateDerivatyX(double x, double y) {
                throw new UnsupportedOperationException("Not supported yet.");
            }

            public double calculateDerivatyY(double x, double y) {
                throw new UnsupportedOperationException("Not supported yet.");
            }
        };


        ArrayList<ArrayList<Particle>> ParticleSupport = generateParticleSupport();
        for (int i = 0; i < ParticleSupport.size(); i++) {
            ListOfShapeFunction.get(i).setNodeList(ParticleSupport.get(i));
        }

        A = new HashMatrix2(HashMatrix2.ARRAY_TYPE, NodeList.size());
        b = new double[NodeList.size()];
        ElementMatrixGeneratorMF2DSQR gen = new ElementMatrixGeneratorMF2DSQR();
        gen.setLeft(left);
        gen.setRight(right);
        gen.setBottom(bottom);
        gen.setTop(top);

        for (int i = 0; i < ParticleSupport.size(); i++) {
//            ArrayList<Particle> arrayList = ParticleSupport.get(i);
            gen.setPhi1(ListOfShapeFunction.get(i));
            for (int j = 0; j < ParticleSupport.size(); j++) {
                gen.setPhi2(ListOfShapeFunction.get(j));
                gen.update();
                double val = gen.generateKeElement(p) + gen.generateMeElement(q) + gen.generatePeElement(w);
                A.setValue(i, j, val);
            }
                b[i]= gen.generateFeElement(f);
        }

        for (int i = 0; i < b.length; i++) {
            for (int j = 0; j < b.length; j++) {
                System.out.print(A.getValue(i, j)+"\t\t\t");
            }
//            System.out.println("\t\t\t\t||\t"+b[i]);
        }
        for(int i = 0; i < b.length; i++) {
            System.out.println("\t\t\t\t||\t"+b[i]);
        }

        System.out.println();
        System.out.println("=========================================================");

    }

    public double[] solvePenalty(){

        double MaxC = 1000000000000.0;
        GaussQuadrature gau = new GaussQuadrature(100);

        Node nodeA = new Node();
        Node nodeB = new Node();

        if(DirihetAx){
            nodeB.x=left;
            nodeA.x=left;
            for(int i=0;i<b.length;i++){
                final JMeshlessInterpolationFunctionLiu f1 = ListOfShapeFunction.get(i);
//                double y = f1.getYii();
//                double tmpVal = f1.calculate(left, y);
                for(int j=0;j<b.length;j++){
                    final JMeshlessInterpolationFunctionLiu f2 = ListOfShapeFunction.get(j);

                    nodeA.y=Math.max(f1.getYii()-f1.getRho(), f2.getYii()-f2.getRho());
                    nodeB.y=Math.min(f1.getYii()+f1.getRho(), f2.getYii()+f2.getRho());
                    nodeB.y=Math.min(top, nodeB.y);
                    nodeA.y=Math.max(bottom, nodeA.y);

//                    double x = this.left, y = Math.max( f1.getYii()-f1.getRho(), f2.getYii()-f2.getRho()),
//                           x1 = this.left, y1 = Math.min( f1.getYii()+f1.getRho(), f2.getYii()+f2.getRho());
                    double val = gau.GaussOnGamma(nodeA, nodeB, new IFunction2D() {

                        public double calculate(double x, double y) {
                            return f1.calculate(x, y)*f2.calculate(x, y);
                        }

                        public double calculateDerivatyX(double x, double y) {
                            throw new UnsupportedOperationException("Not supported yet.");
                        }

                        public double calculateDerivatyY(double x, double y) {
                            throw new UnsupportedOperationException("Not supported yet.");
                        }
                    });
//                    double val = tmpVal*f2.calculate(left, y);
                    A.setValue(i, j,A.getValue(i, j)+ MaxC*(val));
                }
                b[i]+=MaxC*(AxVal);
            }
        }
        if(DirihetBx){
            nodeB.x=right;
            nodeA.x=right;
            for(int i=0;i<b.length;i++){
                final JMeshlessInterpolationFunctionLiu f1 = ListOfShapeFunction.get(i);
//                double y = f1.getYii();
//                double tmpVal = f1.calculate(right, y);
                for(int j=0;j<b.length;j++){
                    final JMeshlessInterpolationFunctionLiu f2 = ListOfShapeFunction.get(j);
                    nodeA.y=Math.max(f1.getYii()-f1.getRho(), f2.getYii()-f2.getRho());
                    nodeB.y=Math.min(f1.getYii()+f1.getRho(), f2.getYii()+f2.getRho());
                    nodeB.y=Math.min(top, nodeB.y);
                    nodeA.y=Math.max(bottom, nodeA.y);

//                    double x = this.left, y = Math.max( f1.getYii()-f1.getRho(), f2.getYii()-f2.getRho()),
//                           x1 = this.left, y1 = Math.min( f1.getYii()+f1.getRho(), f2.getYii()+f2.getRho());
                    double val = gau.GaussOnGamma(nodeA, nodeB, new IFunction2D() {

                        public double calculate(double x, double y) {
                            return f1.calculate(x, y)*f2.calculate(x, y);
                        }

                        public double calculateDerivatyX(double x, double y) {
                            throw new UnsupportedOperationException("Not supported yet.");
                        }

                        public double calculateDerivatyY(double x, double y) {
                            throw new UnsupportedOperationException("Not supported yet.");
                        }
                    });
//                    double val = tmpVal*f2.calculate(right, y);
                    A.setValue(i, j,A.getValue(i, j)+ MaxC*(val));
                }
                b[i]+=MaxC*(BxVal);
            }
        }
        if(DirihetAy){
            nodeB.y=bottom;
            nodeA.y=bottom;
            for(int i=0;i<b.length;i++){
                final JMeshlessInterpolationFunctionLiu f1 = ListOfShapeFunction.get(i);
//                double x = f1.getXii();
//                double tmpVal = f1.calculate(x, bottom);
                for(int j=0;j<b.length;j++){
                    final JMeshlessInterpolationFunctionLiu f2 = ListOfShapeFunction.get(j);
                    nodeA.x=Math.max(f1.getXii()-f1.getRho(), f2.getXii()-f2.getRho());
                    nodeB.x=Math.min(f1.getXii()+f1.getRho(), f2.getXii()+f2.getRho());
                    nodeB.x=Math.min(right, nodeB.x);
                    nodeA.x=Math.max(left, nodeA.x);

//                    double x = this.left, y = Math.max( f1.getYii()-f1.getRho(), f2.getYii()-f2.getRho()),
//                           x1 = this.left, y1 = Math.min( f1.getYii()+f1.getRho(), f2.getYii()+f2.getRho());
                    double val = gau.GaussOnGamma(nodeA, nodeB, new IFunction2D() {

                        public double calculate(double x, double y) {
                            return f1.calculate(x, y)*f2.calculate(x, y);
                        }

                        public double calculateDerivatyX(double x, double y) {
                            throw new UnsupportedOperationException("Not supported yet.");
                        }

                        public double calculateDerivatyY(double x, double y) {
                            throw new UnsupportedOperationException("Not supported yet.");
                        }
                    });
//                    double val = tmpVal*f2.calculate(x, bottom);
                    A.setValue(i, j,A.getValue(i, j)+ MaxC*(val));
                }
                b[i]+=MaxC*(AyVal);
            }
        }
        if(DirihetBy){
            for(int i=0;i<b.length;i++){
                nodeB.y=top;
                nodeA.y=top;
                final JMeshlessInterpolationFunctionLiu f1 = ListOfShapeFunction.get(i);
//                double x = f1.getXii();
//                double tmpVal = f1.calculate(x, top);
                for(int j=0;j<b.length;j++){
                    final JMeshlessInterpolationFunctionLiu f2 = ListOfShapeFunction.get(j);
//                    double x = Math.max( f1.getXii()-f1.getRho(), f2.getXii()-f2.getRho()), y = this.top,
//                           x1 = Math.min( f1.getXii()+f1.getRho(), f2.getXii()+f2.getRho()), y1 = this.top;
                    nodeA.x=Math.max(f1.getXii()-f1.getRho(), f2.getXii()-f2.getRho());
                    nodeB.x=Math.min(f1.getXii()+f1.getRho(), f2.getXii()+f2.getRho());
                    nodeB.x=Math.min(right, nodeB.x);
                    nodeA.x=Math.max(left, nodeA.x);

                    double val = gau.GaussOnGamma(nodeA, nodeB, new IFunction2D() {

                        public double calculate(double x, double y) {
                            return f1.calculate(x, y)*f2.calculate(x, y);
                        }

                        public double calculateDerivatyX(double x, double y) {
                            throw new UnsupportedOperationException("Not supported yet.");
                        }

                        public double calculateDerivatyY(double x, double y) {
                            throw new UnsupportedOperationException("Not supported yet.");
                        }
                    });
//                    double val = tmpVal*f2.calculate(x, top);
                    A.setValue(i, j,A.getValue(i, j)+ MaxC*(val));
                }
                b[i]+=MaxC*(ByVal);
            }
        }

//        for (int i = 0; i < b.length; i++) {
//            for (int j = 0; j < b.length; j++) {
//                System.out.print(A.getValue(i, j)+"\t\t\t");
//            }
//            System.out.println("\t\t\t\t||\t"+b[i]);
//        }
//
//        System.out.println();
//        System.out.println("=========================================================");

            LUFactorization lu = new LUFactorization(A);
        return lu.calculate(b);
    }

    public double[] solve2(){

        double MaxC = 1000000000000.0;
        GaussQuadrature gau = new GaussQuadrature(20);

        Node nodeA = new Node();
        Node nodeB = new Node();
//        ArrayList<JMeshlessInterpolationFunctionLiu> Aprox = new ArrayList<JMeshlessInterpolationFunctionLiu>();
//            for (int i = 0; i < ListOfShapeFunction.size(); i++) {
//                JMeshlessInterpolationFunctionLiu JMeshlessInterpolationFunctionLiu = ListOfShapeFunction.get(i);
//                if (JMeshlessInterpolationFunctionLiu.getXii()==left) {
//                    Aprox.add(JMeshlessInterpolationFunctionLiu);
//                }
//            }
//            for (int i = 0; i < ListOfShapeFunction.size(); i++) {
//                JMeshlessInterpolationFunctionLiu JMeshlessInterpolationFunctionLiu = ListOfShapeFunction.get(i);
//                if (JMeshlessInterpolationFunctionLiu.getXii()==right) {
//                    Aprox.add(JMeshlessInterpolationFunctionLiu);
//                }
//            }

//            for (int i = 0; i < ListOfShapeFunction.size(); i++) {
//                JMeshlessInterpolationFunctionLiu JMeshlessInterpolationFunctionLiu = ListOfShapeFunction.get(i);
//                if ((JMeshlessInterpolationFunctionLiu.getXii()!=left)&&(JMeshlessInterpolationFunctionLiu.getXii()!=right)&&(JMeshlessInterpolationFunctionLiu.getYii()==top)) {
//                    Aprox.add(JMeshlessInterpolationFunctionLiu);
//                }
//            }

//            for (int i = 0; i < ListOfShapeFunction.size(); i++) {
//                JMeshlessInterpolationFunctionLiu JMeshlessInterpolationFunctionLiu = ListOfShapeFunction.get(i);
//                if ((JMeshlessInterpolationFunctionLiu.getXii()!=left)&&(JMeshlessInterpolationFunctionLiu.getXii()!=right)&&(JMeshlessInterpolationFunctionLiu.getYii()==bottom)) {
//                    Aprox.add(JMeshlessInterpolationFunctionLiu);
//                }
//            }


        if(DirihetAx){

            ArrayList<JMeshlessInterpolationFunctionLiu> Aprox = new ArrayList<JMeshlessInterpolationFunctionLiu>();
            for (int i = 0; i < ListOfShapeFunction.size(); i++) {
                JMeshlessInterpolationFunctionLiu JMeshlessInterpolationFunctionLiu = ListOfShapeFunction.get(i);
                if (JMeshlessInterpolationFunctionLiu.getXii()==left) {
//                    Aprox.add(JMeshlessInterpolationFunctionLiu);
                    A.setValue(i, i, MaxC);
                    b[i]=0.0;
                }
            }

            nodeB.x=left;
            nodeA.x=left;
            for(int i=0;i<Aprox.size();i++){
                final JMeshlessInterpolationFunctionLiu f1 = Aprox.get(i);
//                double y = f1.getYii();
//                double tmpVal = f1.calculate(left, y);
                for(int j=0;j<b.length;j++){
                    final JMeshlessInterpolationFunctionLiu f2 = ListOfShapeFunction.get(j);

                    nodeA.y=Math.max(f1.getYii()-f1.getRho(), f2.getYii()-f2.getRho());
                    nodeB.y=Math.min(f1.getYii()+f1.getRho(), f2.getYii()+f2.getRho());
                    nodeB.y=Math.min(top, nodeB.y);
                    nodeA.y=Math.max(bottom, nodeA.y);

//                    double x = this.left, y = Math.max( f1.getYii()-f1.getRho(), f2.getYii()-f2.getRho()),
//                           x1 = this.left, y1 = Math.min( f1.getYii()+f1.getRho(), f2.getYii()+f2.getRho());
                    double val = gau.GaussOnGamma(nodeA, nodeB, new IFunction2D() {

                        public double calculate(double x, double y) {
                            return f1.calculate(x, y)*f2.calculate(x, y);
                        }

                        public double calculateDerivatyX(double x, double y) {
                            throw new UnsupportedOperationException("Not supported yet.");
                        }

                        public double calculateDerivatyY(double x, double y) {
                            throw new UnsupportedOperationException("Not supported yet.");
                        }
                    });
//                    double val = tmpVal*f2.calculate(left, y);
                    A.setValue(i, j,A.getValue(i, j)+ MaxC*(val));
                }
                b[i]+=MaxC*(AxVal);
            }
        }
        if(DirihetBx){
            ArrayList<JMeshlessInterpolationFunctionLiu> Aprox = new ArrayList<JMeshlessInterpolationFunctionLiu>();
            for (int i = 0; i < ListOfShapeFunction.size(); i++) {
                JMeshlessInterpolationFunctionLiu JMeshlessInterpolationFunctionLiu = ListOfShapeFunction.get(i);
                if (JMeshlessInterpolationFunctionLiu.getXii()==right) {
//                    Aprox.add(JMeshlessInterpolationFunctionLiu);
                    A.setValue(i, i, MaxC);
                    b[i]=0.0;
                }
            }
            nodeB.x=right;
            nodeA.x=right;
            for(int i=0;i<Aprox.size();i++){
                final JMeshlessInterpolationFunctionLiu f1 = Aprox.get(i);
//                double y = f1.getYii();
//                double tmpVal = f1.calculate(right, y);
                for(int j=0;j<b.length;j++){
                    final JMeshlessInterpolationFunctionLiu f2 = ListOfShapeFunction.get(j);
                    nodeA.y=Math.max(f1.getYii()-f1.getRho(), f2.getYii()-f2.getRho());
                    nodeB.y=Math.min(f1.getYii()+f1.getRho(), f2.getYii()+f2.getRho());
                    nodeB.y=Math.min(top, nodeB.y);
                    nodeA.y=Math.max(bottom, nodeA.y);

//                    double x = this.left, y = Math.max( f1.getYii()-f1.getRho(), f2.getYii()-f2.getRho()),
//                           x1 = this.left, y1 = Math.min( f1.getYii()+f1.getRho(), f2.getYii()+f2.getRho());
                    double val = gau.GaussOnGamma(nodeA, nodeB, new IFunction2D() {

                        public double calculate(double x, double y) {
                            return f1.calculate(x, y)*f2.calculate(x, y);
                        }

                        public double calculateDerivatyX(double x, double y) {
                            throw new UnsupportedOperationException("Not supported yet.");
                        }

                        public double calculateDerivatyY(double x, double y) {
                            throw new UnsupportedOperationException("Not supported yet.");
                        }
                    });
//                    double val = tmpVal*f2.calculate(right, y);
                    A.setValue(i, j,A.getValue(i, j)+ MaxC*(val));
                }
                b[i]+=MaxC*(BxVal);
            }
        }
        if(DirihetAy){
            ArrayList<JMeshlessInterpolationFunctionLiu> Aprox = new ArrayList<JMeshlessInterpolationFunctionLiu>();
            for (int i = 0; i < ListOfShapeFunction.size(); i++) {
                JMeshlessInterpolationFunctionLiu JMeshlessInterpolationFunctionLiu = ListOfShapeFunction.get(i);
                if ((JMeshlessInterpolationFunctionLiu.getXii()!=left)&&(JMeshlessInterpolationFunctionLiu.getXii()!=right)&&(JMeshlessInterpolationFunctionLiu.getYii()==top)) {
//                    Aprox.add(JMeshlessInterpolationFunctionLiu);
                    A.setValue(i, i, MaxC);
                    b[i]=0.0;

                }
            }
            nodeB.y=bottom;
            nodeA.y=bottom;
            for(int i=0;i<Aprox.size();i++){
                final JMeshlessInterpolationFunctionLiu f1 = Aprox.get(i);
//                double x = f1.getXii();
//                double tmpVal = f1.calculate(x, bottom);
                for(int j=0;j<b.length;j++){
                    final JMeshlessInterpolationFunctionLiu f2 = ListOfShapeFunction.get(j);
                    nodeA.x=Math.max(f1.getXii()-f1.getRho(), f2.getXii()-f2.getRho());
                    nodeB.x=Math.min(f1.getXii()+f1.getRho(), f2.getXii()+f2.getRho());
                    nodeB.x=Math.min(right, nodeB.x);
                    nodeA.x=Math.max(left, nodeA.x);

//                    double x = this.left, y = Math.max( f1.getYii()-f1.getRho(), f2.getYii()-f2.getRho()),
//                           x1 = this.left, y1 = Math.min( f1.getYii()+f1.getRho(), f2.getYii()+f2.getRho());
                    double val = gau.GaussOnGamma(nodeA, nodeB, new IFunction2D() {

                        public double calculate(double x, double y) {
                            return f1.calculate(x, y)*f2.calculate(x, y);
                        }

                        public double calculateDerivatyX(double x, double y) {
                            throw new UnsupportedOperationException("Not supported yet.");
                        }

                        public double calculateDerivatyY(double x, double y) {
                            throw new UnsupportedOperationException("Not supported yet.");
                        }
                    });
//                    double val = tmpVal*f2.calculate(x, bottom);
                    A.setValue(i, j,A.getValue(i, j)+ MaxC*(val));
                }
                b[i]+=MaxC*(AyVal);
            }
        }
        if(DirihetBy){
            ArrayList<JMeshlessInterpolationFunctionLiu> Aprox = new ArrayList<JMeshlessInterpolationFunctionLiu>();
            for (int i = 0; i < ListOfShapeFunction.size(); i++) {
                JMeshlessInterpolationFunctionLiu JMeshlessInterpolationFunctionLiu = ListOfShapeFunction.get(i);
                if ((JMeshlessInterpolationFunctionLiu.getXii()!=left)&&(JMeshlessInterpolationFunctionLiu.getXii()!=right)&&(JMeshlessInterpolationFunctionLiu.getYii()==bottom)) {
//                    Aprox.add(JMeshlessInterpolationFunctionLiu);
                    A.setValue(i, i, MaxC);
                    b[i]=0.0;
                }
            }
                nodeB.y=top;
                nodeA.y=top;
            for(int i=0;i<Aprox.size();i++){
                final JMeshlessInterpolationFunctionLiu f1 = Aprox.get(i);
                
                
//                double x = f1.getXii();
//                double tmpVal = f1.calculate(x, top);
                for(int j=0;j<b.length;j++){
                    final JMeshlessInterpolationFunctionLiu f2 = ListOfShapeFunction.get(j);
//                    double x = Math.max( f1.getXii()-f1.getRho(), f2.getXii()-f2.getRho()), y = this.top,
//                           x1 = Math.min( f1.getXii()+f1.getRho(), f2.getXii()+f2.getRho()), y1 = this.top;
                    nodeA.x=Math.max(f1.getXii()-f1.getRho(), f2.getXii()-f2.getRho());
                    nodeB.x=Math.min(f1.getXii()+f1.getRho(), f2.getXii()+f2.getRho());
                    nodeB.x=Math.min(right, nodeB.x);
                    nodeA.x=Math.max(left, nodeA.x);

                    double val = gau.GaussOnGamma(nodeA, nodeB, new IFunction2D() {

                        public double calculate(double x, double y) {
                            return f1.calculate(x, y)*f2.calculate(x, y);
                        }

                        public double calculateDerivatyX(double x, double y) {
                            throw new UnsupportedOperationException("Not supported yet.");
                        }

                        public double calculateDerivatyY(double x, double y) {
                            throw new UnsupportedOperationException("Not supported yet.");
                        }
                    });
//                    double val = tmpVal*f2.calculate(x, top);
                    A.setValue(i, j,A.getValue(i, j)+ MaxC*(val));
                }
                b[i]+=MaxC*(ByVal);
            }
        }

        for (int i = 0; i < b.length; i++) {
//            for (int j = 0; j < b.length; j++) {
//                System.out.print(A.getValue(i, j)+"\t\t\t");
//            }
            System.out.println("\t\t\t\t||\t"+b[i]);
        }
//
//        System.out.println();
//        System.out.println("=========================================================");

            LUFactorization lu = new LUFactorization(A);
            
//            double[] Ljambda = new double[b.length];
//            for (int i = 0; i < Ljambda.length; i++) {
//                Ljambda[i]=0;
//            }
//            
//            double[] prev = new double[A.getN()];
//            double[] next= null;
//            do {
//            	double[] newB = new double[A.getN()];
//            	int i = 0;
//            	for (Iterator<JMeshlessInterpolationFunctionLiu> iterator = ListOfShapeFunction.iterator(); iterator
//						.hasNext();) {
//            		final JMeshlessInterpolationFunctionLiu funcV = iterator.next();
//            		//upper bound
//            		newB[i] = b[i] + gau.GaussOnGamma(new Node(-0.5, 0.5, true, 0), new Node(0.5, 0.5, true, 1), funcV) + 
//            		gau.GaussOnGamma(new Node(-0.5, -0.5, true, 0), new Node(0.5, -0.5, true, 1), funcV);
//				}
//            	
//            	
//            	next = lu.calculate(newB); 
//			} while (check(prev, next));
            
            
        return lu.calculate(b);
    }
    
    private boolean check(double[] prev, double[] next){
        for(int i=0;i<prev.length;i++){
            if(Math.abs(prev[i]-next[i])>eps){
                System.out.print("перевірка:\t");
                System.out.println(i+"\t"+Math.abs(prev[i]-next[i]));
                return false;
            }
        }
        return true;
    }

    public double[] solveLagrangeMulti(){

        int n = A.getN();

        int num = 0;
        int DirihletAxNum=0;
        int DirihletBxNum=0;
        int DirihletAyNum=0;
        int DirihletByNum=0;
        if(DirihetAx){
            num++;
            DirihletAxNum=n;
        }
        if(DirihetBx){
            DirihletBxNum=n+num*NNN;
            num++;
        }
        if(DirihetAy){
            DirihletAyNum=n+num*NNN;
            num++;
        }
        if(DirihetBy){
            DirihletByNum=n+num*NNN-2;
            num++;
        }

        HashMatrix A2=new HashMatrix(HashMatrix.ARRAY_TYPE, n+num*NNN-4);
        double[] tmpb= new double[n+num*NNN-4];
//перезапис старого ...
            for (int i = 0; i < b.length; i++) {
                for (int j = 0; j < b.length; j++) {
                    A2.setValue(i, j, A.getValue(i, j));
                }
                tmpb[i] = b[i];
            }

        GaussQuadrature gau = new GaussQuadrature(100);
        if(DirihetAx){
            Node nodeA = new Node();
            Node nodeB = new Node();
            nodeA.x=left;
            nodeB.x=left;

            ArrayList<JMeshlessInterpolationFunctionLiu> LagrangeMultiAprox = new ArrayList<JMeshlessInterpolationFunctionLiu>();
            for (int i = 0; i < ListOfShapeFunction.size(); i++) {
                JMeshlessInterpolationFunctionLiu JMeshlessInterpolationFunctionLiu = ListOfShapeFunction.get(i);
                if (JMeshlessInterpolationFunctionLiu.getXii()==left) {
                    LagrangeMultiAprox.add(JMeshlessInterpolationFunctionLiu);
                }
            }


            for(int i=0;i<LagrangeMultiAprox.size();i++){
                final JMeshlessInterpolationFunctionLiu f1 = LagrangeMultiAprox.get(i);
                for(int j=0;j<b.length;j++){
                    final JMeshlessInterpolationFunctionLiu f2 = ListOfShapeFunction.get(j);
                    nodeA.y=Math.max(f1.getYii()-f1.getRho(), f2.getYii()-f2.getRho());
                    nodeB.y=Math.min(f1.getYii()+f1.getRho(), f2.getYii()+f2.getRho());
                    nodeB.y=Math.min(top, nodeB.y);
                    nodeA.y=Math.max(bottom, nodeA.y);
                    if (nodeA.y<nodeB.y) {
//                        double val = gau.GaussOnGamma(nodeA, nodeB, new IFunction2D() {
//
//                            public double calculate(double x, double y) {
//                                return f1.calculate(x, y)*f2.calculate(x, y);
//                            }
//
//                            public double calculateDerivatyX(double x, double y) {
//                                throw new UnsupportedOperationException("Not supported yet.");
//                            }
//
//                            public double calculateDerivatyY(double x, double y) {
//                                throw new UnsupportedOperationException("Not supported yet.");
//                            }
//                        });
                        double val = f2.calculate(f1.getXii(), f1.getYii());
//                        System.out.print((i+n)+"\t"+ j+"\t"+val);
                        A2.setValue(n+i, j, val);
                        A2.setValue(j, n+i, val);
//                        System.out.println("\t\t"+A2.getValue(n+i, j));
                        tmpb[n+i]=gau.GaussOnGamma(nodeA, nodeB, new IFunction2D() {

                            public double calculate(double x, double y) {
                                return f1.calculate(x, y)*AxVal;
                            }

                            public double calculateDerivatyX(double x, double y) {
                                throw new UnsupportedOperationException("Not supported yet.");
                            }

                            public double calculateDerivatyY(double x, double y) {
                                throw new UnsupportedOperationException("Not supported yet.");
                            }
                        });
                    }
                }
            }
        }
        if(DirihetBx){
            Node nodeA = new Node();
            Node nodeB = new Node();
            nodeA.x=right;
            nodeB.x=right;

            ArrayList<JMeshlessInterpolationFunctionLiu> LagrangeMultiAprox = new ArrayList<JMeshlessInterpolationFunctionLiu>();
            for (int i = 0; i < ListOfShapeFunction.size(); i++) {
                JMeshlessInterpolationFunctionLiu JMeshlessInterpolationFunctionLiu = ListOfShapeFunction.get(i);
                if (JMeshlessInterpolationFunctionLiu.getXii()==right) {
                    LagrangeMultiAprox.add(JMeshlessInterpolationFunctionLiu);
                }
            }


            for(int i=0;i<LagrangeMultiAprox.size();i++){
                final JMeshlessInterpolationFunctionLiu f1 = LagrangeMultiAprox.get(i);
                for(int j=0;j<b.length;j++){
                    final JMeshlessInterpolationFunctionLiu f2 = ListOfShapeFunction.get(j);
                    nodeA.y=Math.max(f1.getYii()-f1.getRho(), f2.getYii()-f2.getRho());
                    nodeB.y=Math.min(f1.getYii()+f1.getRho(), f2.getYii()+f2.getRho());
                    nodeB.y=Math.min(top, nodeB.y);
                    nodeA.y=Math.max(bottom, nodeA.y);
                    if (nodeA.y<nodeB.y) {
//                        double val = gau.GaussOnGamma(nodeA, nodeB, new IFunction2D() {
//
//                            public double calculate(double x, double y) {
//                                return f1.calculate(x, y)*f2.calculate(x, y);
//                            }
//
//                            public double calculateDerivatyX(double x, double y) {
//                                throw new UnsupportedOperationException("Not supported yet.");
//                            }
//
//                            public double calculateDerivatyY(double x, double y) {
//                                throw new UnsupportedOperationException("Not supported yet.");
//                            }
//                        });
                        double val = f2.calculate(f1.getXii(), f1.getYii());
//                        System.out.print((i+DirihletBxNum)+"\t"+ j+"\t"+val);
                        A2.setValue(i+DirihletBxNum, j, val);
                        A2.setValue(j, i+DirihletBxNum, val);
//                        System.out.println("\t\t"+A2.getValue(DirihletBxNum+i, j));
                        tmpb[i+DirihletBxNum]=gau.GaussOnGamma(nodeA, nodeB, new IFunction2D() {

                            public double calculate(double x, double y) {
                                return f1.calculate(x, y)*BxVal;
                            }

                            public double calculateDerivatyX(double x, double y) {
                                throw new UnsupportedOperationException("Not supported yet.");
                            }

                            public double calculateDerivatyY(double x, double y) {
                                throw new UnsupportedOperationException("Not supported yet.");
                            }
                        });
                    }
                }
            }
        }
        if(DirihetAy){
            Node nodeA = new Node();
            Node nodeB = new Node();
            nodeA.y=bottom;
            nodeB.y=bottom;

            ArrayList<JMeshlessInterpolationFunctionLiu> LagrangeMultiAprox = new ArrayList<JMeshlessInterpolationFunctionLiu>();
            for (int i = 0; i < ListOfShapeFunction.size(); i++) {
                JMeshlessInterpolationFunctionLiu JMeshlessInterpolationFunctionLiu = ListOfShapeFunction.get(i);
                if ((JMeshlessInterpolationFunctionLiu.getXii()!=left)&&(JMeshlessInterpolationFunctionLiu.getXii()!=right)&&(JMeshlessInterpolationFunctionLiu.getYii()==bottom)) {
                    LagrangeMultiAprox.add(JMeshlessInterpolationFunctionLiu);
                }
            }


            for(int i=0;i<LagrangeMultiAprox.size();i++){
                final JMeshlessInterpolationFunctionLiu f1 = LagrangeMultiAprox.get(i);
                for(int j=0;j<b.length;j++){
                    final JMeshlessInterpolationFunctionLiu f2 = ListOfShapeFunction.get(j);
                    nodeA.x=Math.max(f1.getXii()-f1.getRho(), f2.getXii()-f2.getRho());
                    nodeB.x=Math.min(f1.getXii()+f1.getRho(), f2.getXii()+f2.getRho());
                    nodeB.x=Math.min(right, nodeB.x);
                    nodeA.x=Math.max(left, nodeA.x);
                    if (nodeA.x<nodeB.x) {
//                        double val = gau.GaussOnGamma(nodeA, nodeB, new IFunction2D() {
//
//                            public double calculate(double x, double y) {
//                                return f1.calculate(x, y)*f2.calculate(x, y);
//                            }
//
//                            public double calculateDerivatyX(double x, double y) {
//                                throw new UnsupportedOperationException("Not supported yet.");
//                            }
//
//                            public double calculateDerivatyY(double x, double y) {
//                                throw new UnsupportedOperationException("Not supported yet.");
//                            }
//                        });
                        double val = f2.calculate(f1.getXii(), f1.getYii());
//                        System.out.print((i+DirihletAyNum)+"\t"+ j+"\t"+val);
                        A2.setValue(i+DirihletAyNum, j, val);
                        A2.setValue(j, i+DirihletAyNum, val);
//                        System.out.println("\t\t"+A2.getValue(DirihletAyNum+i, j));
                        tmpb[i+DirihletAyNum]=gau.GaussOnGamma(nodeA, nodeB, new IFunction2D() {

                            public double calculate(double x, double y) {
                                return f1.calculate(x, y)*AyVal;
                            }

                            public double calculateDerivatyX(double x, double y) {
                                throw new UnsupportedOperationException("Not supported yet.");
                            }

                            public double calculateDerivatyY(double x, double y) {
                                throw new UnsupportedOperationException("Not supported yet.");
                            }
                        });
                    }
                }
            }
        }
        if(DirihetBy){
            Node nodeA = new Node();
            Node nodeB = new Node();
            nodeA.y=top;
            nodeB.y=top;

            ArrayList<JMeshlessInterpolationFunctionLiu> LagrangeMultiAprox = new ArrayList<JMeshlessInterpolationFunctionLiu>();
            for (int i = 0; i < ListOfShapeFunction.size(); i++) {
                JMeshlessInterpolationFunctionLiu JMeshlessInterpolationFunctionLiu = ListOfShapeFunction.get(i);
                if ((JMeshlessInterpolationFunctionLiu.getXii()!=left)&&(JMeshlessInterpolationFunctionLiu.getXii()!=right)&&(JMeshlessInterpolationFunctionLiu.getYii()==top)) {
                    LagrangeMultiAprox.add(JMeshlessInterpolationFunctionLiu);
                }
            }


            for(int i=0;i<LagrangeMultiAprox.size();i++){
                final JMeshlessInterpolationFunctionLiu f1 = LagrangeMultiAprox.get(i);
                for(int j=0;j<b.length;j++){
                    final JMeshlessInterpolationFunctionLiu f2 = ListOfShapeFunction.get(j);
                    nodeA.x=Math.max(f1.getXii()-f1.getRho(), f2.getXii()-f2.getRho());
                    nodeB.x=Math.min(f1.getXii()+f1.getRho(), f2.getXii()+f2.getRho());
                    nodeB.x=Math.min(right, nodeB.x);
                    nodeA.x=Math.max(left, nodeA.x);
                    if (nodeA.x<nodeB.x) {
//                        double val = gau.GaussOnGamma(nodeA, nodeB, new IFunction2D() {
//
//                            public double calculate(double x, double y) {
//                                return f1.calculate(x, y)*f2.calculate(x, y);
//                            }
//
//                            public double calculateDerivatyX(double x, double y) {
//                                throw new UnsupportedOperationException("Not supported yet.");
//                            }
//
//                            public double calculateDerivatyY(double x, double y) {
//                                throw new UnsupportedOperationException("Not supported yet.");
//                            }
//                        });
                        double val = f2.calculate(f1.getXii(), f1.getYii());
//                        System.out.print((i+DirihletByNum)+"\t"+ j+"\t"+val);
                        A2.setValue(i+DirihletByNum, j, val);
                        A2.setValue(j, i+DirihletByNum, val);
//                        System.out.println("\t\t"+A2.getValue(DirihletByNum+i, j));
                        tmpb[i+DirihletByNum]=gau.GaussOnGamma(nodeA, nodeB, new IFunction2D() {

                            public double calculate(double x, double y) {
                                return f1.calculate(x, y)*ByVal;
                            }

                            public double calculateDerivatyX(double x, double y) {
                                throw new UnsupportedOperationException("Not supported yet.");
                            }

                            public double calculateDerivatyY(double x, double y) {
                                throw new UnsupportedOperationException("Not supported yet.");
                            }
                        });
                    }
                }
            }
        }


        for (int i = 0; i < tmpb.length; i++) {
            for (int j = 0; j < tmpb.length; j++) {
                System.out.print(A2.getValue(i, j)+"\t\t\t");
            }
            System.out.println("\t\t\t\t||\t"+tmpb[i]);
        }

        System.out.println();
        System.out.println("=========================================================");

        //LUFactorization lu = new LUFactorization(A);
        Gauss lu = new Gauss(A2);
        return lu.calculate(tmpb);
    }

    public double[] solve(){
        HashMatrix A2 = new HashMatrix(HashMatrix.ARRAY_TYPE, A.getN()+32);
        double[] tmpb = new double[A.getN()+32];
        //перезапис старого ...
            for (int i = 0; i < b.length; i++) {
                for (int j = 0; j < b.length; j++) {
                    A2.setValue(i, j, A.getValue(i, j));
                }
                tmpb[i] = b[i];
            }
        int pos = 0;
            for (int i = 0; i < ListOfShapeFunction.size(); i++) {
                JMeshlessInterpolationFunctionLiu JMeshlessInterpolationFunctionLiu = ListOfShapeFunction.get(i);
                if (JMeshlessInterpolationFunctionLiu.getYii()==top) {
                    //A.setValue(i, i, A.getValue(i, i)*100000000000.0);
                    //b[i] = 0.0;
                    for(int j=0;j<b.length;j++){
                    final JMeshlessInterpolationFunctionLiu f2 = ListOfShapeFunction.get(j);
                        double val = f2.calculate(JMeshlessInterpolationFunctionLiu.getXii(), JMeshlessInterpolationFunctionLiu.getYii());
                        A2.setValue(A.getN()+pos, j, val);
                        A2.setValue(j, A.getN()+pos, val);
                        tmpb[A.getN()+pos]=0.0;
                    }
                    pos++;
                    continue;
                }
                if (JMeshlessInterpolationFunctionLiu.getYii()==bottom) {
                   //A.setValue(i, i, A.getValue(i, i)*100000000000.0);
                    //b[i] = 0.0;
                    for(int j=0;j<b.length;j++){
                    final JMeshlessInterpolationFunctionLiu f2 = ListOfShapeFunction.get(j);
                        double val = f2.calculate(JMeshlessInterpolationFunctionLiu.getXii(), JMeshlessInterpolationFunctionLiu.getYii());
                        A2.setValue(A.getN()+pos, j, val);
                        A2.setValue(j, A.getN()+pos, val);
                        tmpb[A.getN()+pos]=0.0;
                    }
                    pos++;
                    continue;
                }
                if (JMeshlessInterpolationFunctionLiu.getXii()==left) {
                    //A.setValue(i, i, A.getValue(i, i)*100000000000.0);
                    //b[i] = 0.0;
                    for(int j=0;j<b.length;j++){
                    final JMeshlessInterpolationFunctionLiu f2 = ListOfShapeFunction.get(j);
                        double val = f2.calculate(JMeshlessInterpolationFunctionLiu.getXii(), JMeshlessInterpolationFunctionLiu.getYii());
                        A2.setValue(A.getN()+pos, j, val);
                        A2.setValue(j, A.getN()+pos, val);
                        tmpb[A.getN()+pos]=0.0;
                    }
                    pos++;
                    continue;
                }
                if (JMeshlessInterpolationFunctionLiu.getXii()==right) {
                    //A.setValue(i, i, A.getValue(i, i)*100000000000.0);
                    //b[i] = 0.0;
                    for(int j=0;j<b.length;j++){
                    final JMeshlessInterpolationFunctionLiu f2 = ListOfShapeFunction.get(j);
                        double val = f2.calculate(JMeshlessInterpolationFunctionLiu.getXii(), JMeshlessInterpolationFunctionLiu.getYii());
                        A2.setValue(A.getN()+pos, j, val);
                        A2.setValue(j, A.getN()+pos, val);
                        tmpb[A.getN()+pos]=0.0;
                    }
                    pos++;
                    continue;
                }
            }


            for (int i = 0; i < tmpb.length; i++) {
            for (int j = 0; j < tmpb.length; j++) {
                System.out.print(A2.getValue(i, j)+"\t\t\t");
            }
            System.out.println("\t\t\t\t||\t"+tmpb[i]);
        }

        System.out.println();
        System.out.println("=========================================================");
        
        Gauss g = new Gauss(A2);
        return g.calculate(tmpb);
    }

    public double getAxVal() {
        return AxVal;
    }

    public void setAxVal(double AxVal) {
        this.AxVal = AxVal;
    }

    public double getAyVal() {
        return AyVal;
    }

    public void setAyVal(double AyVal) {
        this.AyVal = AyVal;
    }

    public double getBxVal() {
        return BxVal;
    }

    public void setBxVal(double BxVal) {
        this.BxVal = BxVal;
    }

    public double getByVal() {
        return ByVal;
    }

    public void setByVal(double ByVal) {
        this.ByVal = ByVal;
    }

    public boolean isDirihetAx() {
        return DirihetAx;
    }

    public void setDirihetAx(boolean DirihetAx) {
        this.DirihetAx = DirihetAx;
    }

    public boolean isDirihetAy() {
        return DirihetAy;
    }

    public void setDirihetAy(boolean DirihetAy) {
        this.DirihetAy = DirihetAy;
    }

    public boolean isDirihetBx() {
        return DirihetBx;
    }

    public void setDirihetBx(boolean DirihetBx) {
        this.DirihetBx = DirihetBx;
    }

    public boolean isDirihetBy() {
        return DirihetBy;
    }

    public void setDirihetBy(boolean DirihetBy) {
        this.DirihetBy = DirihetBy;
    }

    public double getBottom() {
        return bottom;
    }

    public void setBottom(double bottom) {
        this.bottom = bottom;
    }

    public double getTop() {
        return top;
    }

    public void setTop(double top) {
        this.top = top;
    }

    public String getfString() {
        return fString;
    }

    public void setfString(String fString) {
        this.fString = fString;
    }

    public String getpString() {
        return pString;
    }

    public void setpString(String pString) {
        this.pString = pString;
    }

    public String getqString() {
        return qString;
    }

    public void setqString(String qString) {
        this.qString = qString;
    }

    public double getRight() {
        return right;
    }

    public void setRight(double right) {
        this.right = right;
    }

    public String getwString() {
        return wString;
    }

    public void setwString(String wString) {
        this.wString = wString;
    }

    public double getLeft() {
        return left;
    }

    public void setLeft(double left) {
        this.left = left;
    }


    public ArrayList<JMeshlessInterpolationFunctionLiu> getListOfShapeFunction() {
        return ListOfShapeFunction;
    }

    public void setListOfShapeFunction(ArrayList<JMeshlessInterpolationFunctionLiu> ListOfShapeFunction) {
        this.ListOfShapeFunction = ListOfShapeFunction;
    }

    public ArrayList<Particle> getNodeList() {
        return NodeList;
    }

    public void setNodeList(ArrayList<Particle> NodeList) {
        this.NodeList = NodeList;
    }

    private ArrayList<ArrayList<Particle>> generateParticleSupport(){
        ArrayList<ArrayList<Particle>> res = new ArrayList<ArrayList<Particle>>();
        //for (int i = 0; i < ListOfShapeFunction.size(); i++) {
        //    ArrayList<Particle> tmp = new ArrayList<Particle>();
        //    JMeshlessInterpolationFunctionLiu1D phi = ListOfShapeFunction.get(i);
            /*for (int j = 0; j < i; j++) {
            JMeshlessInterpolationFunctionLiu1D phiTMP = ListOfShapeFunction.get(j);
            if(Math.abs(phi.getXii()-phiTMP.getXii())<(Math.abs(phi.getRho()+phiTMP.getRho()))){
            tmp.add(new Particle(phiTMP.getXii(), 0.0, j));//TODO перевірити чи той номер вузла ні на шо не впливає .....
            }
            }*/
        //    for (int j = 0; j < ListOfShapeFunction.size(); j++) {
        //        JMeshlessInterpolationFunctionLiu1D phiTMP = ListOfShapeFunction.get(j);
        //        if(Math.abs(phi.getXii()-phiTMP.getXii())<(Math.abs(phi.getRho()+phiTMP.getRho()))){
        //            tmp.add(new Particle(phiTMP.getXii(), 0.0, j));//TODO перевірити чи той номер вузла ні на шо не впливає .....
        //        }
        //    }
        //    res.add(tmp);
        //}
        //return res;
        for (int i = 0; i < ListOfShapeFunction.size(); i++) {
            res.add(NodeList);
        }
        return res;
    }

    public static void main(String[] args){
        ArrayList<Particle> nodes = new ArrayList<Particle>();
        ArrayList<JMeshlessInterpolationFunctionLiu> listBaseFunction = new ArrayList<JMeshlessInterpolationFunctionLiu>();
        MF2DSolver solver  = new MF2DSolver();
        //параметри області:
        double ax = 0, bx = 1,
               ay = 0, by = 1;

        int n = 7, m = 7;
        solver.NNN = 8;

        double  hx = (bx-ax)/(double)n,
                hy = (by-ay)/(double)m;

        //вузли+функції:
        int iii = 0;
        for (int i = 0; i <= n; i++) {
            for (int j = 0; j <= m; j++) {
                Particle p = new Particle(ax+i*hx, ay+j*hy, iii);
                iii++;
                nodes.add(p);
                listBaseFunction.add(new JMeshlessInterpolationFunctionLiu(1, p.getX(), p.getY(), 1.5*hx, new JGaussian(1.5*hx)));
            }
        }

        
//        граничні умови
        solver.setDirihetAx(true);
        solver.setAxVal(0.0);
        solver.setDirihetBx(true);
        solver.setBxVal(0.0);
        solver.setDirihetAy(true);
        solver.setAyVal(0.0);
        solver.setDirihetBy(true);
        solver.setByVal(0.0);
        //параметри області
        solver.setLeft(ax);
        solver.setRight(bx);
        solver.setBottom(ay);
        solver.setTop(by);
        //параметри рівняння
        solver.setpString("1");
        solver.setqString("0");
        solver.setwString("0");
        solver.setfString("10");
        //базові ф-ції і вузли
        solver.setListOfShapeFunction(listBaseFunction);
        solver.setNodeList(nodes);

        solver.generate();
//        double[] res = solver.solvePenalty();
//        double[] res = solver.solveLagrangeMulti();
//        double[] res = solver.solve();
        double[] res = solver.solve2();

        System.out.println("======================");
        for (int i = 0; i < res.length; i++) {
//            double d = res[i];
            System.out.println(res[i]);
        }
        System.out.println("======================");

        n = 10;
        m = 10;

        hx = (bx-ax)/(double)n;
        hy = (by-ay)/(double)m;

        //вузли+функції:
        for (int i = 0; i <= n; i++) {

            for (int j = 0; j <= m; j++) {
                double tmp = 0.0;
                for (int k = 0; k < listBaseFunction.size(); k++) {
                    JMeshlessInterpolationFunctionLiu JMeshlessInterpolationFunctionLiu = listBaseFunction.get(k);
                    tmp+=res[k]*JMeshlessInterpolationFunctionLiu.calculate(ax+i*hx, ay+j*hy);
                }
                System.out.println((ax+i*hx)+"\t"+(ay+j*hy)+"\t"+tmp);
                //Particle p = new Particle(ax+i*hx, ay+j*hy, i*n+j);
                //nodes.add(p);
                //listBaseFunction.add(new JMeshlessInterpolationFunctionLiu(1, p.getX(), p.getY(), 2, new JGaussian()));
            }
        }

    }
}
